Nonlinear radial oscillations of neutron stars 
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The effects of nonlinear oscillations in compact stars are attracting considerable current interest. 
In order to study such phenomena in the framework of fully nonlinear general relativity, highly 
accurate numerical studies are required. A numerical scheme specifically tailored for such a study 
is based on formulating the time evolution in terms of deviations from a stationary equilibrium 
configuration. Using this technique, we investigate over a wide range of amplitudes nonlinear effects 
in the evolution of radial oscillations of neutron stars. In particular, we discuss mode coupling due 
to nonlinear interaction, the occurrence of resonance phenomena, shock formation near the stellar 
surface as well as the capacity of nonlinearities to stabilize perturbatively unstable neutron star 
models. 

PACS numbers: 04.25.D-, 04.40.Dg, 97.10.Sj 



I. INTRODUCTION 

The observational and theoretical study of stellar oscil- 
lations aimed at gaining insight into stellar structure, an 
area called asteroseismology, has been a rich research field 
for a long time. This includes the period-luminosity re- 
lation for Cepheids and the variability of RR Lyrae stars 
first discovered by Leavitt and Pickering and Flemming, 
respectively, 0, see [E 0, H[ and references therein 
for more recent studies. Similarly, much insight into the 
structure of the Sun has been obtained via helio seismol- 
ogy (see @ for a recent review) . Oscillations have also at- 
tracted a great deal of attention in the context of compact 
objects, neutron stars and black holes; consider, for ex- 
ample, the stability analysis of neutron stars using radial 
oscillation modes by Chandrasekhar [2|. More recently, 
interest in neutron star and black-hole oscillations has 
focused on their potential in the context of gravitational 
wave (GW) physics [8(- 

While most insight into stellar oscillations has been 
obtained from linear analysis (see, for example, [E EE 
ITT1 . [HI), nonlinear effects are known to play an impor- 
tant role in the phenomenology of oscillations in various 
scenarios. Examples include the nonlinear coupling of 
modes in the beat Cepheids (see, for example, [13| ) and 
the saturation of r-mode oscillations and, thus, GW gen- 
eration in rotating neutron stars first studied by Schenk 
et at. [3 and extended in [H EE EE EE El- More 
generally, perturbative studies have been extended to in- 
clude nonlinear effects up to cubic order terms in the 
perturbations 0, IH HH- Mode coupling in neutron 
stars with particular regard to the generation of GWs 
has been investigated in the framework of higher-order 
perturbation theory in [IE HI; see also [25l. IM I27L |2S| 
for a gauge-invariant framework to compute higher-order 
perturbations including fluid backgrounds. While the 
present study focuses on neutron star oscillations, we em- 
phasize that black-hole oscillations play an equally im- 



portant role in GW physics [29|, yfj ■ At the linearized 
level, the correspondingsolutions are given in the form 
of quasinormal modes [H HE HE HH which dominate the 
signal from the late stages of a binary black-hole coales- 
cence (see, for example, [H, HE HE Hzj] ) • Nonlinearities 
in black-hole oscillations have been investigated in the 
framework of higher-order perturbation theory as well 
as numerical relativity [34], Hl| . A conclusive answer re- 
garding the presence of nonlinear signatures in ring-down 
waveforms resulting from binary-black-hole inspiral and 
merger has as yet not been obtained, however, because 
of the high numerical accuracy required for such studies 

In comparison with black-hole configurations gravity 
is about an order of magnitude weaker in spacctimcs 
containing neutron stars. This has motivated a variety 
of hydrodynamic simulations which implement gravita- 
tional effects in the form of some approximation such as 
Newtonian theory, the Cowling approximation or confor- 
mal flatness [H, EE El HE EE El] • The present decade 
has seen dramatic progress in the numerical solution of 
the full Einstein field equations, however, and has re- 
sulted in fully relativistic simulations of single compact 
stars, c ollap se to black holes and neutron star binaries 
EE EE S3 EE HE HE 111 HI HI- As in the case of 
black-hole simulations the available accuracy has not yet 
reached a level to facilitate high-precision studies of non- 
linear effects, in particular in the mildly nonlinear regime. 

In consequence, there currently exists a gap in the lit- 
erature studying in the fully nonlinear general relativis- 
tic framework with high precision mildly and moderately 
nonlinear effects in oscillations of compact stars. The 
main purpose of the present paper is to fill this gap in 
the case of the simplest type of stellar pulsations, ra- 
dial oscillations of spherically symmetric polytropic stel- 
lar models. We are aware that more realistic simulations 
of neutron stars require the inclusion of a multitude of 
microphysical effects such as more realistic equations of 
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state and magnetic fields. For this reason, our study is 
intended in the first place to provide a general taxon- 
omy of nonlinear features that will be encountered in the 
context of stellar pulsations. The results also provide 
valuable tools for calibrating the accuracy of general rel- 
ativistic three-dimensional hydro codes. 

Following Refs. [H, we achieve the necessary ac- 
curacy by formulating the problem in terms of deviations 
from an equilibrium background model, as is commonly 
done in traditional perturbation theory. In contrast to 
that approximation, however, we keep all terms of higher 
order in the deviations and thus arrive at a system of 
equations equivalent to the original nonlinear system. 
The main advantage of this approach is the elimination 
from the equations of all terms exclusively containing 
background quantities and, thus, the discretization er- 
ror associated with these terms. The key improvement 
over the toy problem studied in [H3| is the inclusion of 
the entire star including the important outer layers near 
the stellar surface. 

This paper is organized in five sections. We set up 
the numerical framework in Sec. |TTJ In Sec. IIHI we study 
in detail the coupling of eigenmodes due to nonlinear ef- 
fects, including a more detailed discussion of saturation 
and resonance effects. Section HVl investigates two further 
nonlinear effects not directly concerning mode coupling, 
the stabilization of linearly unstable stars and the forma- 
tion of discontinuities near the stellar surface. Finally, we 
summarize our findings in Sec. |Vj Throughout this work 
we adopt units corresponding to c = G = 1, where c is 
the speed of light and G is the gravitational constant. 



II. NUMERICAL FRAMEWORK 

A. Evolution equations 

The starting point for our description of a dynamic, 
spherically symmetric neutron star is the formulation de- 
veloped by May and White [55|, [56| . Specifically, we write 
the metric in the form 



ds 2 



Ydt 1 + p 2 dx z +r 2 (d6 2 +>. 



(1) 



where A, p and areal radius r are functions of time t and 
radial position x. Unless stated otherwise, the coordinate 
x is initialized by the areal radius of the background con- 
figuration (cf. Sec. lIIBI below) and serves as a Lagrangian 
coordinate following the motion of the fluid elements dur- 
ing the time evolution. The neutron star is modeled as a 
single component perfect fluid at zero temperature. The 
corresponding energy momentum tensor can be expressed 
in the form T» v = (p+P)u^u v +Pg^ v , where p and P are 
energy density and pressure. The four velocity obeys 
the normalization condition u^u^ = — 1. Because the ra- 
dial coordinate x is comoving with the fluid elements we 
have the simple expression = [A -1 , 0,0,0]. We fur- 
ther assume rest mass conservation and neglect all heat 



transfer other than that due to the motion of the fluid. 
In particular, this excludes heat transfer by neutrinos or 
radiation as well as pair production and interaction with 
external fields. The equation of state (EOS) is modeled 
by a polytropic law P = Kp 1 , where K is the polytropic 
constant and 7 the polytropic index. 

In order to achieve high accuracy near the boundaries, 
we find it important to employ variables with at most 
linear asymptotic behavior at both the center and the 
stellar surface 1 . We therefore describe the stellar model 
in terms of a rescaled mass 2 function N and a function 
a which are defined by 



N 



1 f l 1 
2r 1 



P 

a = — 
P 



(2) 
(3) 



The set of variables is completed by the metric compo- 
nent A and an auxiliary velocity function w = r tt /X intro- 
duced to obtain a system of partial differential equations 
of first order in space and time. 

The Einstein equations G Mt , = 87rTp„ and the conser- 
vation of energy and momentum V^T^ = then result 
in the following equations: 
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In these equations, commas are used to represent partial 
derivatives and C 2 = dp/ dp is the sound speed. Only 
five of these equations are independent and we choose to 
use Eq. (J9j) in the determination of the eigenmodes below 
but not in the nonlinear time evolution. 

In order to determine boundary conditions, we first 
consider the origin. Because the physical system is spher- 
ically symmetric by construction, all functions that are 
odd in the radius necessarily vanish at the origin, that 
is r(0) = 0, N(0) = and w(0) = 0. The surface a; s 
of the star is defined as the radius where the pressure 
vanishes. We further define our time coordinate such 



1 See Sec. V B.l in |57|| for a discussion of the difficulties arising 
from higher-order falloff of variables near the boundaries. 

2 A straightforward calculation shows that m = r 2 N evaluated 
at the surface of the star is the mass parameter of the exterior 
vacuum metric in Schwarzschild coordinates. 
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that the interior solution matches to the Schwarzschild 
metric in the exterior. We thus obtain outer boundary 
conditions P(xs) = 0, corresponding to cr(irg) = 0, and 
X(x s ) = y/l-2N(x s )r(x s ). 



B. Deviations from an equilibrium configuration 

The purpose of this study is to analyze nonlinear effects 
in the time evolution of radial oscillations of a neutron 
star model with particular regard to the mildly nonlinear 
regime. Such an investigation requires numerical simu- 
lations of high accuracy. We achieve this by formulating 
the numerical evolution in terms of deviations from a 
stationary background configuration. In such a formula- 
tion all terms involving exclusively background quantities 
and, thus, the discretization errors associated with these 
terms drop out of the equations. Because we are mostly 
interested in scenarios where the deviations are small in 
comparison with the background terms, the elimination 
of these error terms leads to a significant increase in nu- 
merical accuracy. 

This decomposition requires us to choose a convenient 
background. For the study at hand, the obvious choice is 
a spherically symmetric, static neutron star, as described 
by Tolman, Oppenheimer and Volkoff (TOV) [H, HI]. 
The system of equations governing the TOV equilibrium 
configuration is given by the time independent limit of 
the system of Eqs. (g])-©, 



A 



N x 

X jX 



X 

r, x ( 47rp- 2 — 



A 



1 - 2Nf 



(N + AnPr), 



(10) 
(11) 
(12) 



where we have introduced an overbar to distinguish back- 
ground variables from their time dependent counterparts. 
Note that f x = 1 if we initialize x by the areal radius 
of the background model. Next, we formally introduce 
deviations from this equilibrium by 



f(x,t)=f(x)+Af(x,t). 



(13) 



Here and in the following, the function / stands for either 
of the variables {A, N, a, p, P, r}. Following customary 
notation in the literature, we denote the radial displace- 
ment by £ = Ar. 

We can now insert Eq. (fl3|) into the set of evolution 
Eqs. (H])-©. Of particular significance in the resulting 
expressions are the terms containing exclusively back- 
ground quantities and no deviations. A straightforward 
calculation demonstrates that all these terms drop out of 
the equations because, by construction, the background 
variables / obey the stationary limit of Eqs. (HJ)-©, that 
is, the TOV equations. We are left exclusively with terms 



of linear or higher order in the deviations A/ and the 
evolution equations are given by 

a. x {AXC 2 + AAC* 2 ) AA, X (1 + a){C 2 - a) 



Act . 



AC* 2 AC 2 
\ x [Aa{C 2 - a) + (1 + a)(AC 2 - Act)] 
AC* 2 



AN, X = f x 4wAp + ^ x 4np- 

N,xt + Cx2N + f x (2AN - 4ttp0 
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2ANr -2N^) — 



— (NX + iirPrX) - — (ANX - AXN) 



^4tt [P(f AA + £A) + APrA] 



(18) 



The boundary conditions for the deviations follow di- 
rectly from those obtained for the full variables / in 
Sec. MM They are given by £(0) = 0, AiV(0) = 0, 
w(0) = as well as Aa(x s ) = and AA(A + AA) = 
-2(£iV + fAN) . 



C. Eigenmodes 

In the limit of infinitcsimally small deviations Af from 
the equilibrium solution /, the time evolution of the star 
is governed by the linearized version of the system (fT4|) - 
(1181) . The solutions to the linearized perturbation equa- 
tions are given by an infinite set of oscillation modes with 
characteristic frequencies. This set of modes is a vital as- 
set when interpreting the nonlinear time evolution. For 
each neutron star model, we therefore calculate the first 
10 eigenmode profiles and their associated frequencies. 

The linearized evolution system is most conveniently 
written in terms of the rescaled radial displacement £ = 
f 2 £/A (cf. chapter 26 in [6(|) and reduces to the single 
partial differential equation 



QC 



(19) 



Here, the variables W, H and Q depend only on the back- 
ground quantities. Exact expressions for these auxiliary 
functions are given in Eqs. (|A1[) - (|A3[) . 

Separation of variables straightforwardly implies a har- 
monic time dependence, i. e. C ~ e lult , where w is a fre- 
quency yet to be determined. The resulting equation 
represents a singular Sturm-Liouville problem, and the 
solutions C^^) °f this eigenvalue problem form a com- 
plete orthonormal set. It is possible, therefore, to expand 
the radial displacement function Q(t, x) in a series 



C(t,x) = J2MtM° 



(20) 
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where the eigenmode coefficients Ai(t) are given by the 
inner product 



A i (t) = (C,Ci)= WC(t,x)d(x) dx . (21) 
Jo 

The Ai(t) thus represent a measure of how strongly an 
eigenmode i contributes at a given time t. 



D. Numerical methods 

The equations for the TOV background ((T0 |i - fpL2 |) as 
well as the Sturm-Liouville problem (|19p represent two- 
point boundary-value problems which can be solved by 
means of a relaxation method [6l| . We discretize the evo- 
lution system l|14[) - (fTg)) using the implicit, second-order 
accurate Crank-Nicholson scheme. From a numerical 
point of view, this implicit evolution scheme turns out 
to be conceptually identical to the two-point boundary- 
value problems of the background and eigenmode calcu- 
lation. We therefore use the same relaxation algorithm 
for all of these calculations. 

Finally, we compute the integrals appearing in the cal- 
culation of the inner products using a standard fourth- 
order accurate Simpson's Rule algorithm. 



E. Code tests 

The code has been tested in several ways. First, we 
have performed a convergence analysis for a polytropic 
model with 7 = 2 and K = 150 km 2 . We use 801, 1601 
and 3201 grid points and obtain second-order conver- 
gence for the background model, the calculation of the 
eigenmode profiles and the time evolution including the 
eigenmode coefficients. In Fig. [TJ we show the conver- 
gence factors of the £2 norms of the evolution variables 
Act, AN and £ for a simulation of a perturbation given 
initially in the form of the first eigenmode with an am- 



plitude of 10 m. The convergence factor is defined as 
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FIG. 1: The convergence factor q of Act, AN and £ as func- 
tions of time. The obtained results are consistent with second- 
order convergence corresponding to q = 4. 
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(22) 



where the subscripts c, m and / stand for coarse, medium 
and fine resolution. For a second-order scheme and the 
above-mentioned grid setups we expect q — 4. The minor 
deviations from this value are due to zero crossings of the 
variables involved. 

For the second test, we use the constraint Eq. (|9|) 
which vanishes in the continuum limit. In order to quan- 
tify the constraint violation due to the discretization er- 
ror, we normalize the numerical constraint by the sum of 
the £2 norms of the individual terms. We find the result- 
ing normalized constraint violations to be of the order of 
10~ 6 or less. 

Third, we test the eigenmode calculation by check- 
ing the completeness of the eigenmode spectrum. This 
is done by calculating the weighted sum X^^iACiO) 
which has to be unity due to the completeness of the ba- 
sis. In our simulations we include the first 10 terms in this 
infinite series and obtain a maximal deviation from the 
expected value of 1 by about 10~ 4 . As a by-product, this 
result demonstrates that the first 10 eigenmodes capture 
most of the dynamics of the system. 

Finally, we compare the eigenmode frequencies with 
results available in the literature and observe excellent 
agreement; the maximal deviation in the first three eigen- 
modes is about 1% from results reported in [l^] and less 
than 1% from those of [62| . 

Using the setup described in this section, we obtain 
highly accurate time evolutions. In fact, we observe that 
the overall error is dominated by the calculation of the 
inner products in Eq. (|2ip via the Simpson algorithm, 
rather than the time evolution of the grid variables. We 
believe this to be an artifact of the uncertainties in the 
eigenmode profiles themselves, the relative error being 
of the order of 10 -5 . Motivated by this assumption, 
we managed to further improve the accuracy using the 
fact that our time evolutions are typically dominated by 
one particular eigenmode. First, we identify this mode 
from the initial data construction. In the course of the 
evolution, we calculate the eigenmode coefficient associ- 
ated with this mode. Before calculating the other eigen- 
mode coefficients, however, we subtract the contribution 
of the dominating mode from the grid variables. By 
virtue of the orthogonality of the eigenmodes, this sub- 
traction does not affect the other coefficients in the con- 
tinuum limit. For finite numerical resolutions, however, 
the eigenmodes are not perfectly orthogonal and we avoid 
contamination of the overlap integrals due to the domi- 
nant mode. 

Taking into account all numerical effects, we arrive at 
the following estimates for the uncertainties. All simula- 
tions discussed in the next section start from initial data 
consisting of a single mode. This mode is found to domi- 
nate the ensuing evolution and we measure its amplitude 
with a relative error of 10~ 5 . All other modes are ab- 
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sent in the initial data, but are excited due to nonlinear 
effects. The accuracy of their measurement depends on 
their amplitude. For those modes used in our analysis we 
obtain relative uncertainties ranging from 10~ 2 for weak 
excitation to 10~ 5 for strong excitation. 



ILL MODE COUPLING 

Before analyzing in detail the coupling of eigenmodes 
in our simulations, we need to address a conceptual diffi- 
culty arising from the nonlinear nature of hydrodynam- 
ics. Eigenmodes are, by construction, a feature of a linear 
theory and our calculation of the eigenmode spectrum in 
Sec. Ill CI required us to specify a background configu- 
ration. Given a dynamic system evolved with the fully 
nonlinear theory, we analyze the deviations from that 
background configuration by projecting them onto the 
eigenmodes associated with that same background. The 
problem is that there exists no unique way of decom- 
posing such a dynamic system into a background plus 
perturbations. Different background configurations im- 
ply different eigenmode spectra. In consequence of this 
ambiguity, we can interprete nonzero mode coefficients 
for modes not present in the initial data in two ways; (i) 
excitation due to nonlinear effects and (ii) a change in the 
background configuration and its eigenmode spectrum 
due to the finite amplitude of the initial perturbation. 
Both interpretations reflect the nonlinear nature of the 
theory and the nonuniqueness of choosing a background 
makes it impossible to distinguish in a well-defined man- 
ner between those interpretations. In the remainder of 
this work we will use the terminology mode coupling and 
excitation but, as a reminder of the ambiguities in the 
interpretation, we will put the words in italics. Note, 
however, that once a particular choice of the background 
has been specified, the ensuing projection of deviations 
from that background onto the corresponding eigenmode 
spectrum is uniquely defined. In that sense, the following 
analysis is self-consistent and any study agreeing upon 
the same decomposition can be compared straightfor- 
wardly with our results. 

Specifically, our analysis is based on the following con- 
struction of a background model plus nonlinear devia- 
tions. We consider either of the two stellar models labeled 
1 and 2 in Table [T] Numerically, we prescribe the initial 
perturbation of this background model by displacing the 
fluid elements from their equilibrium positions according 
to the profile of the displacement £ corresponding to one 
single eigenmode. From this displacement, we straight- 
forwardly calculate the profiles for the remaining evolu- 
tion variables AA, Act and AN according to Eqs. (fT5)) . 
(H3J) and (|15[) in that order. Finally, we set the radial 
velocity w to zero. In the remainder of this work, we 
always label the single mode present in the initial data 
as j and denote its amplitude in the form of its surface 
displacement £j = £(xs) Similarly, we use the label i to 
identify all modes present in the time evolution. 



Model 


7 


K 


Pc [g/cm 3 ] 


M [M ] 


R [km] 


1 


2 


150 km 2 


2.022 x 10 15 


1.555 


10.82 


2 


2.25 


700 km 2 5 


3.600 x 10 15 


1.584 


8.414 


3 


2 


150 km 2 


3.774 x 10 15 


1.655 


9.222 



TABLE I: Physical parameters for the three stellar models 
studied in this work. 

During the evolution, we measure the presence of a 
given mode i in the form of the associated eigenmode 
coefficient Ai(t). In practice, we find these coefficients 
to oscillate periodically; cf. Figures [5] and [3] It is conve- 
nient, therefore, to use the resulting oscillation maximum 
Ai = max. \Ai(t)\ as an overall measure of the degree of 
excitation of that mode in the entire evolution. The in- 
tuitive interpretation of the resulting maxima is greatly 
helped by their one-to-one relation with the correspond- 
ing displacement of the stellar surface. We recall for this 
purpose Eq. (|20[) . By setting Ak = 1 and Ai = for all 
i =/= k in that relation, we can evaluate the surface dis- 
placement C( x s) and, thus, £(xs) corresponding to eigen- 
mode k with unit amplitude Ak- In Table HT1 we list the 
resulting surface displacements for both stellar models 
studied in this section. It becomes clear from the table 
that we expect to deal with numerically small values of 
the coefficients Ak <C 1. 

In the remainder of this section we will analyze in detail 
how initial data given in the form of one single eigenmode 
gives rise to the excitation of other modes. 

A. Model 1: A 7 = 2 polytrope 

1. Exciting the fundamental mode 

We first focus on the coupling of eigenmodes in the case 
of model 1, a 7 = 2 polytrope which is stable against (lin- 
ear) radial perturbations. Our first analysis is based on 
an initial displacement of the fluid elements given by the 
profile of the fundamental eigenmode j = 1. The simula- 
tions last for a physical time of 5 ms and the amplitude 
of the initial data is varied between 10 cm and 70 m. 
For illustration, we show in Fig. [5] the time evolution of 
the eigenmode coefficients A\ , . . . , A4 obtained for an ini- 
tial surface displacement £1 = 10 m. Higher-order modes 
show a similar behavior with decreasing amplitude. Note 
that the coefficients Ai corresponding to weakly excited 
modes such as i = 3 and i = 4 in the figure do not os- 
cillate around zero but reveal an offset. This effect is 
not surprising bearing in mind that the sinusoidal de- 
pendence of the eigenmode oscillations is a direct conse- 
quence of the separation of variables in Eq. (|19[) accord- 
ing to C{t> x ) — Ci( x ) elUit - 111 the linear limit we expect 
an eigenmode to describe a stellar oscillation around the 
equilibrium configuration. In the present case, however, 
the weak eigenmodes have to be viewed as perturbations 
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FIG. 2: The evolution of the first four eigenmode coefficients A% obtained for model 1 and an initial displacement corresponding 
to the fundamental mode with a surface amplitude of £1 = 10 m 



of a time dependent background, namely the equilibrium 
configuration plus any strongly excited eigenmodes. The 
observed offset thus arises as one of the nonlinear effects 
of the system. Indeed, we empirically find that it in- 
creases quadratically with the amplitude of the initial 
data and disappears as we reduce the initial displacement 
of the fluid elements to zero. 

The excitation of higher-order modes inevitably corre- 
sponds to a flow of energy away from the initially present 
mode. Whereas this effect is too small in the example of 
Fig.[Uto be noticeable in the amplitude A\(t), the strong 
coupling between modes 2 and 4 in the case of model 2 
in Fig. [3] demonstrates a strong amplitude modulation of 




time [ms] 



FIG. 3: The envelopes of the amplitudes of the second and 
fourth eigenmode coefficients Ai obtained for model 2 and an 
initial displacement corresponding to the second mode with a 
surface amplitude of £2 = 10 m. Note that the excitation of 
mode i — 4 is so strong that it leads to a visible temporary 
decrease in A2 fu- 



tile mode coefficients. We will discuss this exceptionally 
strong coupling in more detail in the context of resonance 
effects in Sec. MI Dl below. 

As mentioned above, we use the maxima for each of 
the |Aj(t)| to determine the degree of excitation of the 
individual modes. The resulting values for the first six 
eigenmodes are shown in Fig. [4] We have ignored higher- 
order modes because they are very weakly excited and 
therefore subject to uncertainties larger than the limit 
discussed in Sec. Ill El 

The results in the figure show that the amplitude of 
the fundamental eigenmode grows linearly with the initial 
amplitude £1, as expected. In contrast, all other eigen- 
mode coefficients exhibit a quadratic dependence on £1 

At = . (23) 

The expansion coefficients c^i obtained from regression 
are listed in Table HTT1 Throughout the entire amplitude 
range we do not observe any signs of a transition of the 
power laws from quadratic to higher order. This indicates 
that the nonlinear interaction is dominated by leading- 
order effects. This changes when we consider different 
initial data. 



2. Exciting higher modes 

Next, we consider initial perturbations in the form of 
the second and the third eigenmode, respectively. The 
corresponding plots, now up to and including mode 10, 
are shown in Fig. [5] As before the eigenmode coefficient 
of the initially present mode shows linear behavior. For 
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Mode 



£(a;s) [km] model 1 



£(xs) [km] model 2 



Mode 



£(xs) [km] model 1 



£(xs) [km] model 2 



2.15 
6.31 
11.08 
16.41 
22.27 



I. 54 
3.84 
6.29 
8.88 

II. 60 



6 
7 
8 
9 
10 



28.64 
35.47 
42.74 
50.42 
58.49 



14.45 
17.43 
20.51 
23.70 
26.98 



TABLE II: Surface displacement £(x s ) corresponding to eigenmode k with unit amplitude Ak = 1 for the first two stellar models 
of Table U 




ft [km] 

FIG. 4: The maximal eigenmode coefficients of the first six 
eigenmodes as functions of the amplitude of the initial per- 
turbation in form of the fundamental mode. The curves rep- 
resent a linear function in £i for the fundamental mode and 
quadratic power laws for all other modes. 



i 2 3 4 5 6 

Ci,i[km~ 2 }4.9 x 10~ 3 2.9 x 10~ a 7.9 x 10~ 4 1.3 x 10~ 4 1.1 x 10~ 5 



TABLE III: The quadratic expansion coefficients for the mode 
coupling between the fundamental and other modes. 

sufficiently small amplitudes of the initial displacement, 
we also observe the quadratic power law for all modes 
not present in the initial data. For larger perturbations, 
however, these show a clear transition to higher-order 
power laws. This transition demonstrates a significant 
contribution of couplings beyond the leading-order terms 
and we will refer to this regime as the moderately non- 
linear regime as opposed to the weakly nonlinear regime 
where quadratic coupling terms dominate. The details of 
the transition such as the threshold amplitude of the ini- 
tial data and the relative significance of the higher-order 
terms depend on the individual mode under considera- 
tion as well as on the choice of initial data. We empiri- 
cally investigate the transition by fitting the eigenmode 
coefficients according to 

^ = c i j&) 2 +d i j&) 3 + e id fa)* + fiAZi)*+-'- ■ ( 24 ) 

In practice, we find this series to be dominated by a sub- 
set of the terms on the right-hand side and we explicitly 



set subdominant coefficients to zero. The resulting coef- 
ficients are listed in Table IIVI and the corresponding fits 
are shown as the lines in Figs. [4] [5] and [6j While there 
appears to be a tendency for the expansion coefficients 
Ci.j, di.j, ... to decrease for larger mode numbers i, there 
are exceptions to this rule. Consider, for example, C3.2 
and C4 ; 2 which demonstrate that the second mode cou- 
ples more strongly to mode 4 than to mode 3. This is 
corroborated by the data for mode 3 (filled diamonds) 
and 4 (filled upward triangle) in the left panel of Fig. [5j 
It is interesting to note in this context that the frequency 
ratio of modes 2 and 4 is close to 1:2. We will return to 
the issue of resonance phenomena in Sec. IIII Dl 



Another feature visible both in the expansion coeffi- 
cients in Table IIVI and the curves in Fig. [5] is that the 
modes excited due to nonlinear effects form groups of 
similar polynomial behavior. For j = 2 they form pairs 
and for j = 3 triplets. As a function of the initial am- 
plitude of mode 2 (mode 3), for example, modes 5 and 6 
(modes 7, 8 and 9) show a transition to a cubic power law 
and modes 7 and 8 (modes 10, 11 and 12) a transition 
to a fourth-order power law. See also the block diagonal 
structure in Table ITVl listing the coefficients di$, and 
fi,2 (coefficients di^ and e^). 

To summarize our findings, we observe a weakly non- 
linear regime in which the amplitude of secondary modes 
grows quadratically with the amplitude of the initial dis- 
placement. In the case of j = 2 or j = 3 and sufficiently 
large amplitudes, the mode excitation exhibits a tran- 
sition to higher-order power laws. Secondarily excited 
modes form multiplets; for initial data in the form of 
mode j > 1, the excitation of modes i = j + l,...,2j 
depends quadratically on the initial amplitude, modes 
i — 2j + 1, 3j show a transition to a third-order power 
law, modes i = 3j + l, 4j a transition to a fourth-order 
power law and so on. Finally, higher-order modes appear 
to have a significantly stronger tendency to transfer en- 
ergy to lower-order modes than the other way around. 
The latter observation was also indicated by Papadopou- 
los and Sopuerta's [38] study of black-hole oscillations 
and suggested to explain the robustness of black-holes to 
strong deformations. 
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FIG. 5: The maximal eigenmode coefficients of the first ten eigenmodes as functions of the amplitude of the initial perturbation 
in form of the second (left panel) and the third (right panel) mode. 



i 


Ci,2 




di,2 


e;,2 


fi.2 


c;,3 




di,3 


ei,3 




[km" 2 ] 




[km" 3 ] 


[km" 4 ] 


[km" 5 ] 


[km- 2 ] 




[km" 3 ] 


[km" 4 ] 


1 


0.63 













0.64 










2 












0.095 










3 


0.033 





















4 


0.12 













0.020 










5 


4.1 x 10" 


-4 


0.078 








0.029 










6 


1.3 x 10" 


-4 


0.069 








0.029 










7 


1.9 x 10" 


-5 





0.118 





5.8 x 10" 


4 


0.040 





8 


9.4 x 10" 


-6 





0.073 





4.9 x 10" 


5 


0.034 





9 


9.3 x 10" 


-6 








0.172 


2.3 x 10" 


5 


0.019 





10 


1.1 x 10" 


-5 








0.101 


1.1 x 10" 


5 





0.057 


11 












1.1 x 10" 


5 





0.036 


12 












1.2 x 10" 


5 





0.020 



TABLE IV: Expansion coefficients for model 1 obtained from fitting Eq. (I24|l to the data points for initial data in the form of 
mode j = 2 and j = 3, respectively. 



B. Model 2: A stiffer equation of state 



Saturation 



We now turn our attention to model 2 of Table Q] We 
want to assess to what extent the observations in the pre- 
vious section depend on the equation of state of the stel- 
lar model. As before, we consider three scenarios where 
we prescribe initial data in the form of mode j = 1, 2 or 3 
and measure the degree of excitation of the other modes 
up to and including i — 10 due to nonlinear effects. It is 
reassuring to see that using this model we qualitatively 
reproduce most of the general features discussed above. 
For comparison with the results for model 1, we show in 
Tables |V] and IVII the expansion coefficients and in Fig. [6] 
the eigenmode coefficients resulting from initial data in 
the form of eigenmode j = 2 or 3. 

While confirming the overall behavior observed for 
model 1, however, the analysis of model 2 reveals addi- 
tional complications arising from two effects which were 
barely visible in the analysis of model 1: saturation and 
resonance. 



We first discuss the saturation effect observed in the 
analysis of model 2. Consider for this purpose the left 
panel of Fig. [6] which shows the excitation of modes re- 
sulting for j — 2. We first notice a relatively strong 
excitation of mode 4; compare the filled upper triangles 
in this plot with their counterparts in Fig. [5] for model 
1. Indeed, the excitation of mode 4 is so strong that the 
resulting power law fit crosses the curve for mode 2 at 
initial amplitudes just above lO -2 km. Note, however, 
how the actual mode coefficients for mode 4 (filled up- 
per triangles) start deviating from the polynomial fit at 



i [km' 2 ] 0.032 0.011 1.5 ■ 10~ 3 7.6 ■ 10~ 5 8.9 ■ 10" 



TABLE V: The quadratic expansion coefficients Cj,i for the 
mode coupling between the fundamental and other modes for 
the stiffer neutron star model. 
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FIG. 6: The maximal eigenmode coefficients of the first ten eigenmodes as functions of the amplitude of the initial perturbation 
in form of the second eigenmode (left panel) or the third eigenmode (right panel) as obtained for model 2. 



i 


a, 2 [km 2 ] di,2 


[km" 3 ] 


ei,2[km 4 ] 




d,3 [km -2 ] 


di,3 [km -3 ] 


ei,3 [km 4 ] 


1 


2.8 











3.14 








2 










0.22 








3 


0.09 

















4 


24 











0.053 








5 


2.1 x 10 -3 


21 








0.084 








6 


5.0 x 10 -4 


31 








0.099 








7 


2 x 10~ 6 





2210 





1.3 x 10 -3 


0.22 





8 


< 2 x 10~ 6 





1650 





4.1 x 10~ 4 


0.18 





9 


< 2 x 10" 6 





62.8 





2.0 x 10 -4 


0.11 





10 


< 2 x 10" 6 





35.2 





1.5 x 10 -4 





0.52 



TABLE VI: Expansion coefficients for model 2 obtained from fitting Eq. (|24|l to the data points for initial data in the form of 
mode j = 2 and 3, respectively. 



such large amplitudes; the actual excitation of mode 4 is 
weaker than expected from extrapolation of the polyno- 
mial fit and the coefficients |A4| max remain below those 
of the initially present mode |A2| m ax- It appears that at 
sufficiently large amplitudes, mode 4 starts acting as a 
significant source of excitation itself and transfers energy 
to other modes via nonlinear coupling. The equilibrium 
reached in this way sets a limit for the growth in am- 
plitude of mode 4. On closer investigation of the left 
panel of Fig. [51 we notice a similar behavior for modes 
5 to 8. These modes show saturation at smaller ampli- 
tudes which is compatible with our general observation 
that higher-order modes tend to transfer energy to lower- 
order modes more efficiently than the other way round. 



We conclude from this observation that for the present 
example the assumption of two-mode coupling would 
oversimplify the situation. An accurate modeling of the 
chain of excitation whereby energy is transferred from 
mode 2 via mode 4 to other modes requires us to take 
into account higher-order effects in the mode coupling as 
provided, for example, by our fully nonlinear numerical 
treatment. 



D. Resonance effects 

Resonance is a well-known phenomenon arising in the 
context of the forced oscillator. Following Van Hoolst 
(23 | . we can model the stellar oscillations as a system of 
forced oscillators. The time evolution of the oscillation 
amplitudes is then determined by 

— —t^- +ujfAi — —7 + (a n sin nClt + b n cos nClt) . (25) 
at 2 - ' 2 '-^ 

n 

Here f2 is the frequency and the a n and b n are the ampli- 
tudes of the driving forces. Assuming for simplicity that 
the b n vanish, Eq. ([25| has the analytic solution 

AAt) = V 9 Q " Q sinnttt n e Af . (26) 
' uif — m!r 

n 

In our case, the driving terms are typically dominated 
by the single eigenmode j present in the initial data and 
therefore f2 = Uj. This mode should excite with particu- 
lar efficiency those modes which have an eigen frequency 
close to an integer multiple of u>j . In Table I VIII we show 
the frequencies and corresponding ratios for some modes 
of model 2 whose frequencies have ratios close to integer 
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Model 1 


Model 2 


Model 1 


Model 2 


Mode i 


u)i [kHz] 


u)i [kHz] 






1 


10.908 


11.197 


- 


- 


2 


35.299 


47.526 


1.975 2.890 


2.000 2.939 


3 


53.011 


72.024 


1.925 


1.940 


4 


69.711 


95.065 


1.919 


1.931 


5 


85.985 


117.527 






6 


102.026 


139.694 






7 


117.935 


161.685 






8 


133.756 


183.570 







initial amplitude 
[m] 


frequency fi 
[Hz] 


0.1 


8.21 


1 


8.27 


10 


12.29 


20 


18.98 


30 


25.30 



TABLE VIII: This table presents the frequencies obtained by 
fitting the envelopes of the eigenmode coefficients evolution 
with sinusoidal functions. 



TABLE VII: The frequencies and the most promising reso- 
nance factors L02i/iOi or u^i/iOi for some of the eigenmodes of 
model 1 and model 2. 



values. Resonance occurs most conspicuously between 
modes 2 and 4 and we have already noted the strong ex- 
citation of mode 4 in the left panel of Fig. [5] The effect is 
not as dramatic in the case of modes 6 and 8, but we still 
observe a preferred excitation of these modes compared 
with model 1; compare, for example, their excitation with 
that of mode 3 and the corresponding results in the left 
panel of Fig. [5l We also note a small deviation in the 
transition to different power laws in the moderately non- 
linear regime. Consider modes i — 9, 10 in the case of 
j = 2 in Table IVIl According to the rule, the excitation 
of these modes should exhibit a transition to fifth order 
in £2, but instead we observe a fourth-order dependence. 
We believe this to be a consequence of the strong exci- 
tation of mode 4, so that our approximation of a single 
strong mode driving the coupling is no longer valid. 

In all cases of mode coupling, we observe a periodic 
transfer of energy back and forth between the eigenmodes 
involved. This manifests itself in a modulation of the os- 
cillation amplitude as apparent, for example, in Fig. [2] in 
particular the upper right panel. Such amplitude mod- 
ulation is also present in the case of strong resonance. 
In contrast to "normal" coupling of eigenmodes, how- 




50 100 150 200 250 
time [ms] 



FIG. 7: The evolution of the fourth eigenmode coefficient for 
different perturbation amplitudes. 



ever, we further observe in cases of strong resonance a 
modulation period which depends sensitively on the am- 
plitude of the initial data. We illustrate this in Fig. [7| 
which displays the coefficient A±(t) obtained for model 2 
and initial data in the form of the second eigenmode with 
different amplitudes. As is evident in the figure, larger 
amplitudes result in shorter modulation periods. The re- 
sulting numerical values of the modulation frequency are 
shown in Table I Villi and are well approximated by the 
linear relation 



fimodta) = 7.55 Hz + 577 Hz 



km 



(27) 



In the perturbative limit £2 - > we obtain a modulation 
frequency of fi mo d = 7.55 Hz. A detailed interpreta- 
tion in the context of coupled harmonic oscillators [22j 
is beyond the scope of this paper, largely because it is 
highly nontrivial to calculate the coupling constants. We 
note, however, the close resemblance between the mea- 
sured modulation frequency and the difference between 
the involved eigenmodes' frequencies 



6u> 

2" 



C1J4 — 2 • U>2 



6.5 Hz 



(28) 



Note that 5uj is orders of magnitude below the indi- 
vidual frequencies 0J2 and lo^ and therefore subject to 
a larger relative error. Within this error, the measured 
value agrees with the prediction of Eq. (f2"7|) . 

Leaving a detailed investigation for future work, we 
tentatively interpret our observations as follows. The 
modulation frequency is closely related to the difference 
in the frequencies of the coupling modes and decreases 
significantly as we approach resonance. As has already 
been discussed in Sec. IIII C[ nonlinear effects set a limit 
on the resonance. They further result in a deviation of 
the modulation frequency from the limit of perfect reso- 
nance £! mo d -> Hz. 

In comparison, the resonance between modes 2 and 4 
of model 1 is less pronounced and we find the modula- 
tion frequency to be independent of the amplitude of the 
initial data. Indeed, Eq. (f2"5|) predicts Suj — 141 Hz for 
model 1. This value is an order of magnitude larger than 
that obtained for model 2 and appears to be much more 
robust to effects of a finite initial amplitude £2 as given 
by the second term on the right-hand side of Eq. (I27p . 
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IV. FURTHER NONLINEAR EFFECTS 

In this section we consider two nonlinear effects which 
are not directly related to the coupling of eigenmodes. 
First, we study the stability properties of a stellar model 
close to the maximum of the mass-radius relation but 
located on the unstable branch. From linear theory, we 
would expect this model to be unstable to small pertur- 
bations away from its equilibrium configuration. Second, 
we investigate the consequences of the vanishing of the 
speed of sound at the stellar surface with regard to the 
formation of discontinuities near the surface. 



A. Stabilization of linearly unstable stars 

A particularly interesting scenario for the investigation 
of nonlinear effects is that of marginally stable neutron 
stars. It is well known that neutron stars of sufficient 
compactness become unstable with respect to their fun- 
damental radial oscillation mode [7]. The onset of this 
instability occurs at zero frequency of the fundamental 
mode. At this point linear perturbation terms cancel in 
the equations and the higher-order perturbations become 
more important. 

We study this effect quantitatively for model 3 of Table 
HI We have already mentioned that this model is unsta- 
ble, i.e. its fundamental oscillation mode has an imag- 
inary frequency lo\ = —779 Hz 2 . We test the predic- 
tion of the linearized theory by prescribing initial data in 
the form of this fundamental cigcnmodc using an ampli- 
tude of 10 m. The surprising result of the fully nonlin- 
ear time evolution is shown in Fig. [5J the star oscillates 
periodically over many milliseconds without any sign of 
instability. We find this general behavior independent of 
whether the initial displacement represents an expansion 
(dashed curve) or a compression (solid curve) of the star, 
even though these two scenarios differ in the details of 
the ensuing oscillations. This result demonstrates that 
nonlinear effects have the capacity to stabilize a linearly 
unstable neutron star. It will be interesting in future 
work to study this effect in more detail from an analytic 
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FIG. 8: The coefficient of the fundamental eigenmode. The 
dashed line corresponds to an initial expansion of the star, 
the solid line to a compression. 




FIG. 9: The frequency of the fundamental mode as function of 
the initial perturbation. The frequencies have been obtained 
by a Fourier analysis of the nonlinear evolution. 



point of view. 

In addition to this qualitative difference between per- 
turbative predictions and the nonlinear evolution, we ob- 
serve a variety of quantitative deviations. Contrary to ex- 
pectations from linearized theory, the sign of the initial 
perturbation of the star has an impact on the resulting 
oscillation pattern. In particular, the frequency and am- 
plitude of the oscillation differ significantly for the two 
cases as is apparent in Fig. [H] We further observe an off- 
set of the eigenmode coefficient indicating that the star 
is no longer oscillating symmetrically around its equilib- 
rium position. We have already identified such an offset 
as a nonlinear effect in Sec. IIII A II In the present exam- 
ple, however, the offset is much larger. We further note 
that the solid curve does not have sinusoidal shape but 
is visibly distorted. Finally, we find that the oscillation 
frequency of the nonlinear evolution depends not only on 
the sign but also on the initial amplitude of the pertur- 
bation. Indeed, the measured values differ substantially 
from those predicted by linear theory This is illustrated 
in Fig. [5J where we show the oscillation frequency as a 
function of the amplitude in the range —40 m to +40 m. 
As expected intuitively, we obtain larger real frequen- 
cies, i. e. larger deviations from the linear prediction, for 
larger amplitudes of the initial data. Conversely, we re- 
cover the linear limit and observe collapse of the stellar 
model when choosing initial perturbations of sufficiently 
small amplitude, in the decimeter range for the present 
example. 



B. Shock formation at the surface 

Analysis of the asymptotic structure of the TOV 
model combined with a polytropic equation of state 
(cf. Sec. 5.2.5 of (62j) shows that the speed of sound 
vanishes as (fs — f) 1 / 2 at the surface regardless of the 
parameters of the polytrope. Here fs is the radius of 
the star. Because of this decrease of the speed of sound, 
a signal of finite width propagated towards the surface 
will be compressed; its tail moves faster than its head. 
Naturally one may ask under which circumstances this 
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FIG. 10: Snapshots of the evolution of the velocity w for a moderate amplitude B — 2 x 10 4 of the initial perturbation (left 
panel) and of Act obtained for an amplitude B = 10 -3 . 



gives rise to shock formation. In order to investigate this 
question, we consider model 1 of Table HI 

First, however, we need to take care of a numerical 
difficulty arising in this context. In order to guarantee 
adequate numerical resolution of the expected features 
close to the surface, we switch from the radial coordinate 
x to a rescaled radius y related by 



r tV — C , 



(29) 



where C 2 = is the speed of sound. A straightfor- 
ward calculation shows that the speed of sound measured 
in terms of y now approaches a finite value at the sur- 
face, thus eliminating the danger of under resolving fea- 
tures 3 . Even more remarkably, this decompactification 
is achieved by mapping a finite interval in x to a finite 
interval in the new variable y. For more details on this 
procedure we refer the reader to [H[ and [12] . For the nu- 
merical analysis we perturb model 1 by setting the initial 
velocity 



w{t = 0) = Be-to- V0 ?/ h , 



(30) 



where width and center of the Gaussian pulse are set to 
b = 0.05 km 2 and y — 30 km, respectively. The gen- 
eral behavior of this pulse during a fully nonlinear time 
evolution now depends exclusively on its initial ampli- 
tude B. In practice, we observe three different types of 
behavior. For sufficiently small amplitudes, the pulse is 
reflected at the surface and we do not observe shock for- 
mation. At moderate amplitudes the pulse gets reflected 
at the boundary, but a discontinuity in w as well as Act 
develops shortly after reflection. This is shown in the left 
panel of Fig. [TU] for an amplitude B = 2 x 10~ 4 . As we 
further increase the amplitude, discontinuities develop in 



the variables w, AA and Act before the pulse reaches the 
stellar surface. This scenario is illustrated for the vari- 
able Act and using an initial amplitude B = 10~ 3 in the 
right panel of Fig. [101 

Shock formation at the stellar surface has recently been 
studied in the context of a plane parallel model in Newto- 
nian gravity by Gundlach and Please [64]. In particular, 
their Eq. (27) provides the minimum amplitude required 
for shock formation in terms of amplitude vo as well as 
width ct i and initial location cto of the initial pulse. (We 
use a tilde to distinguish their a from our evolution vari- 
able defined in Eq. ([3]).) This equation neglects a factor 
1/2, which we reintroduce for the comparison with our 
numerical simulations. The starting point for the follow- 
ing discussion is therefore 



vp_ < 1 



2 Uo 



(31) 



where the polytropic index is defined by 7 = 1 + 1/n. 
Translating this expression into the coordinates used 
throughout this work, we arrive at the following condi- 
tion 4 



D 



y/g = const. 



(32) 



Here = |rg — fo\, fo is the position and f\ the charac- 
teristic length scale of the initial perturbation and rs and 
g are radius and gravitational acceleration at the surface 
of the background model, respectively. The derivation 
of relation (f3~lj) does not take into account the precise 
shape of the wave packet, however, which may enter as a 
dimensionless constant. 

In order to test their prediction numerically we perform 
a scries of simulations with Gaussian initial data. As a 



3 We note that this potential lack of numerical resolution is not 
a problem in the evolution of the eigenmodes discussed above 4 This exp ression corrects Eq. (28) in [6J by adding a missing 
because of their oscillatory character. factor *J x+Jxo. 
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+ 3/2 




in [km] 


for initial 


width y 


/ &72 




0.2 


0.3 


0.4 


0.5 


6.79 


0.180 


0.196 


0.204 


0.207 


7.13 


0.168 


0.182 


0.189 


0.191 


7.45 


0.157 


0.170 


0.176 


0.178 


7.76 


0.148 


0.160 


0.165 


0.167 


8.05 


0.141 


0.152 


0.157 


0.158 


8.32 


0.135 


0.145 


0.150 


0.151 


8.58 


0.130 


0.140 


0.144 


0.145 


8.82 


0.126 


0.135 


0.139 


0.140 


9.05 


0.123 


0.132 


0.135 


0.136 


9.27 


0.120 


0.129 


0.132 


0.133 


9.47 


0.118 


0.126 


0.130 


0.130 


9.66 


0.117 


0.125 


0.128 


0.129 


9.83 


0.115 


0.123 


0.127 


0.128 


9.99 


0.114 


0.123 


0.126 


0.127 


10.13 


0.114 


0.122 


0.126 


0.128 


10.27 


0.114 


0.122 


0.126 


0.128 


10.38 


0.114 


0.123 


0.127 


0.130 


10.46 


0.114 


0.123 


0.129 


0.131 


10.58 


(0.115) (0.125) (0.131) (0.134) 


10.65 


(0.116) (0.127) (0.134) (0.139) 


10.72 


(0.118) (0.130) (0.138) (0.141) 



TABLE IX: The left-hand side of Eq. J32J for critical pertur- 
bations that are sufficiently large to generate shocks close to 
the surface. 



measure of the characteristic length scale f\ in Eq. (|32|) . 
we choose the full width at half maximum (FWHM) of the 
Gaussian. For our model 1 we have a polytropic index 
n = 1 and y/g = y/M/r'g = 0.1337. The resulting two- 
parameter study is summarized in Table IIXI For each 
combination of width f\ and initial position f we start 
with a very small amplitude of the initial perturbation 
and repeat the simulation at increasingly larger ampli- 
tude until we observe formation of a discontinuity. The 
resulting threshold amplitude vo is then used to evaluate 
the left-hand side of Eq. ([52]). 

The numerical uncertainties in this study are domi- 
nated by the difficulties in determining when a discon- 
tinuity has actually formed. First, we have to discretize 
the amplitude vq in order to keep the number of simu- 
lations at an acceptable level. Second, numerical dissi- 
pation may suppress shock formation. Finally, there is 
some freedom in choosing a critical gradient to define a 
discontinuity. In our simulations we choose the critical 
value to be w >y = 0.05 km . Bearing in mind all these 
difficulties, we estimate the numerical accuracy of the 
critical amplitude v to be about 10 %. Within these 
uncertainties, we observe good agreement with the pre- 
dictions of [6J]. The obtained values for D are listed in 
Table El 



As long as the initial position of the pulse is not too 
close to the center the left-hand side of Eq. (|32|) shows 
little variation with the pulse width and initial location. 
Bearing in mind that (|32[) has been derived for plane par- 
allel geometry and in the Newtonian approximation, the 
observed increase in variability of the coefficient in Table 
IIXI at small xq is not surprising. For radii close to the sur- 
face the determination of the width and the position of 
the pulse is limited because the pulse steepens and devi- 
ates significantly from a purely Gaussian profile. Values 
obtained in this regime are thus given in parentheses. 

In summary, Eq. (|32[) as well as the numerical study 
illustrate nicely that the significance of nonlinear effects 
does not only depend on the amplitude of the evolv- 
ing feature. Instead, the length scale on which profiles 
change is equally important and may generate signifi- 
cant deviations from perturbative predictions at compar- 
atively small amplitudes of the signal. 

For completeness we note that the stellar model used 
in this study oversimplifies the structure of the stellar 
surface. A crust expected to form at the surface of a 
neutron star is likely to result in more complex phenom- 
ena and requires a more detailed treatment. As has been 
mentioned above, however, the primary purpose of our 
study is to probe the taxonomy of nonlinear effects in 
the weakly and moderately nonlinear regime. We post- 
pone the sophistication of a realistic neutron star model 
to future work. In this spirit, our analysis leads to two 
main conclusions. First, the theoretical modeling of neu- 
tron star surfaces requires particular care and, second, 
the surface exhibits a rich phenomenology, not always 
concurrent with immediately intuitive expectations. 



V. SUMMARY 

In this paper we have investigated a variety of nonlin- 
ear effects associated with radial oscillations of polytropic 
neutron star models. We have performed our study in the 
framework of fully nonlinear general relativity and have 
achieved a high level of numerical precision by formulat- 
ing the time evolution in terms of finite deviations from 
an equilibrium configuration which in our case is given 
by a TOV model. 

The first nonlinear effect studied in detail is the cou- 
pling of eigenmodes due to nonlinear effects. In the ab- 
sence of a unique decomposition of a fully nonlinear dy- 
namic model into background and deviations, we have 
performed one-parameter studies using as initial data the 
equilibrium configuration perturbed in the form of one 
single eigenmode profile with varying amplitude. Our 
results and their interpretation are to be understood in 
the context in this particular construction of a reference 
background. By virtue of the completeness of the eigen- 
mode spectrum, we are able to measure the excitation 
of modes not contained in the initial data by calculat- 
ing eigenmode coefficients from overlap integrals. This 
analysis has revealed two qualitatively different regimes 
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depending on the amplitude of the initial data. For suf- 
ficiently low values, we find the excitation of modes to 
increase quadratically with the amplitude. In the con- 
text of analytic studies based on coupled oscillators, as 
for example employed by Van Hoolst (22J, we interpret 
this quadratic dependence as a leading-order coupling be- 
tween different modes. In the case of an initial perturba- 
tion given by the fundamental eigenmode, this quadratic 
dependence persists over the entire range of initial surface 
displacements of up to 70 m. Initial data of higher modes 
with sufficient amplitude, however, result in a transition 
in the excitation of modes to higher-order power laws. 
Specifically, our results indicate that secondarily excited 
modes appear in multiplcts; initial data of eigenmode j 
excites modes i = j + 1, ...,2j according to a quadratic 
power law, modes i = 2j + 1, ...,3j with cubic depen- 
dence and so on. Strong resonance effects may cause 
some complications, however, which manifest themselves 
in deviations from this rule. This behavior demonstrates 
the significance of higher-order coupling of eigenmodes. 
The onset of this moderately nonlinear regime occurs at 
amplitudes of the order of a few meters, but the details 
appear to vary with the stellar model. As a general rule, 
we find higher-order modes to pass energy to lower-order 
modes more efficiently than the other way around. 

Our results confirm the intuitive expectation that 
eigenmodes with an integer frequency ratio interact par- 
ticularly efficiently. This resonance manifests itself most 
conspicuously in the second and fourth eigenmode of the 
stiffer model considered in this study where the frequency 
ratio is equal to two within three digits. Initial data in the 
form of the second mode can excite the fourth mode with 
such efficiency that both amplitudes become comparable. 
When this happens, we observe a further complication in 
the nonlinear behavior: the strong excitation of mode 4 
results in a significant transfer of energy to other eigen- 
modes and its amplitude saturates, i.e. it stops growing 
in accordance with the expected power law. A similar 
behavior is observed for other pairs of modes, as, for ex- 
ample, modes 2 and 6. Overall, resonance appears to 
be weaker for larger deviations of the involved frequency 
ratios from an integer value. 

A particularly fertile ground for the analysis of nonlin- 
ear effects is given by stellar models close to the stability 
limit of the mass-radius relation. Instability manifests 
itself in the vanishing of the oscillation frequency of the 
fundamental eigenmode 0] in which case the linear (in 
the deviations) terms in the evolution equations cancel. 
The resulting dominance of higher-order terms leads to 
a variety of effects. Most notably, our results show that 
nonlinear effects have the capacity to stabilize stars ex- 
pected to be unstable in linear analysis; initial perturba- 
tions lead to a periodic pulsation of the star instead of 
the expected gravitational collapse. Nonlincarity further 
manifests itself in a visible distortion of the time depen- 
dence of the oscillations away from sinusoidal character. 
Also, the oscillation pattern depends on the initial phase 
of the perturbation. It will be interesting to investigate in 



future work to what extent this stabilization is a generic 
feature and how it depends on the physical properties of 
the stellar model. 

Finally, we have studied the formation of discontinu- 
ities near the surface of the star by evolving Gaussian 
pulses propagating from the stellar interior towards the 
surface. In particular, we have compared our numeri- 
cal findings with analytic predictions by Gundlach and 
Please [64| who give threshold amplitudes for the initial 
amplitude of the pulse depending on its initial location 
in the sense that amplitudes above this threshold lead 
to shock formation and those below do not. Within the 
accuracy of the analysis, approximately 10 %, our nu- 
merical results confirm the predictions of Gundlach and 
Please in the range of validity of their model. 

In summary, the precision of the numerical method 
proposed here has enabled us to identify a number of 
nonlinear effects with no approximation other than that 
of a rather simple stellar model. Our results suggest a 
variety of extensions for future work. These include the 
refinement of the stellar model by including, for exam- 
ple, realistic equations of state, magnetic fields and, most 
importantly, relaxing the condition of spherical symme- 
try. The main obstacle for the latter appears to us to be 
the extension of the comoving, Lagrangian formulation 
to higher dimensional problems. We find a careful treat- 
ment of the stellar surface to be crucial in our study. A 
possible avenue towards achieving this goal is the use of 
level set methods and fast marching methods J65| specif- 
ically designed to follow the motion of interfaces in nu- 
merical simulations. 
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APPENDIX A: THE EIGENMODE EQUATION 

In the limit of infinitcsimally small deviations A/ from 
the equilibrium variables /, the time evolution of the 
star is governed by the linearized version of the system 
(P3J) - (fT5|) . The combination of these equations into the 
single, second-order partial differential Eq. (|19l) for the 
rescaled displacement £ is discussed in detail in [60( • The 
coefficients W , II and Q in that equation are given by 

n = c 2 ( P - + p)^ (Ai) 
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W 



^,x\ 4A X 0-25 



(A2) 
(A3) 



where the metric function ft = (1 — 2Af)~ 1 / 2 . 

In order to numerically solve the cigenmode equation, 
we write it as a first-order system in the variable £ and 
make the ansatz £(t,x) — £(x)e luJt . 



= C 2 i x 



B 



= (uj 2 ), x 



(A4) 



7=5 — X +r x (o; 2 W + Q)k 



Jj* (A5) 
(A6) 




The last equation simply expresses the fact that the 
frequency u> is constant. The purpose of introducing 
this equation becomes apparent when we consider the 
boundary conditions. Spherical symmetry requires the 
displacement £ to vanish at the origin. The auxiliary 
variable B vanishes by definition at the stellar surface. 
The value of B at the origin determines the amplitude of 
the profile £(x). Because eigenmodes are only defined up 
to a constant factor, we are free to choose for B(0) an ar- 
bitrary finite value. We thus have a two-point-boundary- 
value problem for the three variables £(#), B{x) and lo(x) 
which we solve with a relaxation algorithm. In particu- 
lar, the solution for the constant u>(x) provides us with 
the frequency of the eigenmode. Solutions only exist for 
specific values of the frequency and finding a given mode 
requires a relatively careful initial guess. 
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